#1
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
pbmc.markers %>% group_by(cluster) %>% top_n(n = 3, wt = avg_log2FC)

# find markers for every cluster compared to all remaining cells, report only the positive ones
> pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
Calculating cluster 0
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=01s  
Calculating cluster 1
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02s  
Calculating cluster 2
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=01s  
Calculating cluster 3
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=01s  
Calculating cluster 4
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=01s  
Calculating cluster 5
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02s  
Calculating cluster 6
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02s  
Calculating cluster 7
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03s  
Calculating cluster 8
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=01s  


> head(pbmc.markers)
              p_val avg_log2FC pct.1 pct.2     p_val_adj cluster  gene
RPS12 1.273332e-143  0.7298951 1.000 0.991 1.746248e-139       0 RPS12
RPS6  6.817653e-143  0.6870694 1.000 0.995 9.349729e-139       0  RPS6
RPS27 4.661810e-141  0.7281575 0.999 0.992 6.393206e-137       0 RPS27
RPL32 8.158412e-138  0.6196246 0.999 0.995 1.118845e-133       0 RPL32
RPS14 5.177478e-130  0.6252832 1.000 0.994 7.100394e-126       0 RPS14
RPS25 3.244898e-123  0.7496479 0.997 0.975 4.450053e-119       0 RPS25


> pbmc.markers %>% group_by(cluster) %>% top_n(n = 3, wt = avg_log2FC)
# A tibble: 27 × 7
# Groups:   cluster [9]
       p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene  
       <dbl>      <dbl> <dbl> <dbl>     <dbl> <fct>   <chr> 
 1 3.75e-112       1.09 0.912 0.592 5.14e-108 0       LDHB  
 2 9.57e- 88       1.36 0.447 0.108 1.31e- 83 0       CCR7  
 3 1.35e- 51       1.08 0.342 0.103 1.86e- 47 0       LEF1  
 4 0               5.57 0.996 0.215 0         1       S100A9
 5 0               5.48 0.975 0.121 0         1       S100A8
 6 3.89e-268       4.55 1     0.516 5.34e-264 1       LYZ   
 7 1.06e- 86       1.27 0.981 0.643 1.45e- 82 2       LTB   
 8 3.44e- 59       1.22 0.651 0.245 4.71e- 55 2       CD2   
 9 2.97e- 58       1.23 0.42  0.111 4.07e- 54 2       AQP3  
10 0               4.31 0.936 0.041 0         3       CD79A 
# … with 17 more rows
#  Use `print(n = ...)` to see more rows



# fig1
VlnPlot(pbmc, features = c("MS4A1", "CD79A"))

# you can plot raw counts as well
VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)
VlnPlot(pbmc, features = c("NKG7", "PF4"), pt.size = 0)


# fig2：默认使用 slot=data 的数据。
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP",
                               "CD8A"))
#2.2 除了基因，还可以画 PC_1, UMAP_1，及 meta.data 表头: head(pbmc@meta.data)
FeaturePlot(object = pbmc, features = c('PC_1', "UMAP_1", "nFeature_RNA", "percent.mt"), 
            cols = c("#eeeeee", "red") )

#2.3 对于2个基因，还支持blend=T模式
FeaturePlot(pbmc, features = c("CD4", "CD8A"), blend = T)
#2.4 order=T 防遮挡
FeaturePlot(pbmc, features = c("CD4", "CD8A"), blend = T, order = T)

#2.5 使用 split.by 参数分面
pbmc$gene_high=ifelse(pbmc$nFeature_RNA>1500, T, F)
FeaturePlot(object = pbmc, features = c('PC_1', "CD3G"), 
            split.by = "gene_high",
            cols = c("#eeeeee", "red") )
#2.6 标记 ident ，并指定颜色、字号
FeaturePlot(pbmc, features = c("CD4"), label = T, label.size = 5, label.color = "red")


# fig3 heatmap
top5 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 5, wt = avg_log2FC)
DoHeatmap(pbmc, features = top5$gene, label = T)  #+ NoLegend()


DotPlot(pbmc, features = unique(top5$gene)) + RotatedAxis()



# tree
pbmc_tree <- BuildClusterTree(object = pbmc)
Seurat::PlotClusterTree(pbmc_tree)
all.markers <- FindAllMarkers(object = pbmc_tree, node = 14)

top5 <- all.markers %>% group_by(cluster) %>% top_n(n = 5, wt = avg_log2FC)
head(top5)
# fig4
DotPlot(pbmc, features = unique(top5$gene), cluster.idents = T) + RotatedAxis()



- VlnPlot 见 解析8
- FeaturePlot todo
- DoHeatmap todo
- DotPlot todo





()
$ find . | grep -P "(R|c)$"  | xargs grep -in "FindAllMarkers" --color=auto 2>/dev/null
./seurat-4.1.0/R/differential_expression.R:45:FindAllMarkers <- function(


#' Gene expression markers for all identity classes
#'
#' Finds markers (differentially expressed genes) for each of the identity classes in a dataset
#'
#' @inheritParams FindMarkers
#' @param node A node to find markers for and all its children; requires
找到该node及其子节点的标记基因。

#' \code{\link{BuildClusterTree}} to have been run previously; replaces \code{FindAllMarkersNode}
#' @param return.thresh Only return markers that have a p-value < return.thresh, or a power > return.thresh (if the test is ROC)
仅返回P<该值的标记基因。 如果是ROC检验，power>该值才返回该基因。

#'
#' @return Matrix containing a ranked list of putative markers, and associated
#' statistics (p-values, ROC score, etc.)
返回矩阵，包含排序后的推测的标记基因，相关统计量（p值，ROC打分等）

#'
#' @importFrom stats setNames
#'
#' @export
#'
#' @aliases FindAllMarkersNode 别名
#' @concept differential_expression
#'
#' @examples
#' data("pbmc_small")
#' # Find markers for all clusters
#' all.markers <- FindAllMarkers(object = pbmc_small)
#' head(x = all.markers)
#' \dontrun{
#' # Pass a value to node as a replacement for FindAllMarkersNode
#' pbmc_small <- BuildClusterTree(object = pbmc_small)
#' all.markers <- FindAllMarkers(object = pbmc_small, node = 4)
#' head(x = all.markers)
#' }
#'
FindAllMarkers <- function(
  object,
  assay = NULL,
  features = NULL,
  logfc.threshold = 0.25,
  test.use = 'wilcox',
  slot = 'data',
  min.pct = 0.1,
  min.diff.pct = -Inf,
  node = NULL,
  verbose = TRUE,
  only.pos = FALSE,
  max.cells.per.ident = Inf,
  random.seed = 1,
  latent.vars = NULL,
  min.cells.feature = 3,
  min.cells.group = 3,
  pseudocount.use = 1,
  mean.fxn = NULL,
  fc.name = NULL,
  base = 2,
  return.thresh = 1e-2,
  densify = FALSE,
  ...
) {
  #(A1) 自定义函数：功能用到再探究
  MapVals <- function(vec, from, to) {
    vec2 <- setNames(object = to, nm = from)[as.character(x = vec)]
    vec2[is.na(x = vec2)] <- vec[is.na(x = vec2)]
    return(unname(obj = vec2))
  }

  #(A2) 如果是roc检验，默认是wilcox，跳过这里
  if ((test.use == "roc") && (return.thresh == 1e-2)) {
    return.thresh <- 0.7
  }

  #(A3) 如果 node 是空，一般不传入tree，那么就是空
  if (is.null(x = node)) {
  	#(B1) 获取每个细胞的idents属性: 
  	#> length(Idents(pbmc))
	#[1] 2638
	#> head(Idents(pbmc))
	#AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1 AAACCGTGTATGCG-1 AAACGCACTGGTAC-1 
	#               2                3                2                1                6                2 
	#Levels: 0 1 2 3 4 5 6 7 8
	#
	# 取unique值，并排序
	#> sort(x = unique(x = Idents(object = pbmc)))
	#[1] 0 1 2 3 4 5 6 7 8
	#Levels: 0 1 2 3 4 5 6 7 8
    idents.all <- sort(x = unique(x = Idents(object = object)))
  } else {
  	#(B2) 如果传入的是tree，才走这里:跳过
    if (!PackageCheck('ape', error = FALSE)) {
      stop(cluster.ape, call. = FALSE)
    }
    tree <- Tool(object = object, slot = 'BuildClusterTree')
    if (is.null(x = tree)) {
      stop("Please run 'BuildClusterTree' before finding markers on nodes")
    }
    descendants <- DFT(tree = tree, node = node, include.children = TRUE)
    all.children <- sort(x = tree$edge[, 2][!tree$edge[, 2] %in% tree$edge[, 1]])
    descendants <- MapVals(
      vec = descendants,
      from = all.children,
      to = tree$tip.label
    )
    drop.children <- setdiff(x = tree$tip.label, y = descendants)
    keep.children <- setdiff(x = tree$tip.label, y = drop.children)
    orig.nodes <- c(
      node,
      as.numeric(x = setdiff(x = descendants, y = keep.children))
    )
    tree <- ape::drop.tip(phy = tree, tip = drop.children)
    new.nodes <- unique(x = tree$edge[, 1, drop = TRUE])
    idents.all <- (tree$Nnode + 2):max(tree$edge)
  }

  #(A4) 声明两个列表
  genes.de <- list()
  messages <- list()

  #(A5) 遍历每个idents分类(重点)
  for (i in 1:length(x = idents.all)) {
  	#(B1) 默认 verbose=T，开始分析第一个分组
    if (verbose) {
      message("Calculating cluster ", idents.all[i])
    }
    #(B2) 计算 该列表中的第i项
    genes.de[[i]] <- tryCatch(
      expr = {
      	# 调用函数 FindMarkers，前一篇讲过，这里从略
        FindMarkers(
          object = object,
          assay = assay,
          # 如果node为空，则ident.1=ident[i]
          ident.1 = if (is.null(x = node)) {
            idents.all[i]
          } else {
            tree
          },
          # 如果node为空，则不设置 ident.2: 默认就是和其余细胞比。
          ident.2 = if (is.null(x = node)) {
            NULL
          } else {
            idents.all[i]
          },
          features = features,
          logfc.threshold = logfc.threshold,
          test.use = test.use,
          slot = slot,
          min.pct = min.pct,
          min.diff.pct = min.diff.pct,
          verbose = verbose,
          only.pos = only.pos,
          max.cells.per.ident = max.cells.per.ident,
          random.seed = random.seed,
          latent.vars = latent.vars,
          min.cells.feature = min.cells.feature,
          min.cells.group = min.cells.group,
          pseudocount.use = pseudocount.use,
          mean.fxn = mean.fxn,
          fc.name = fc.name,
          base = base,
          densify = densify,
          ...
        )
      },
      #错误：返回错误信息
      error = function(cond) {
        return(cond$message)
      }
    )
    #如果是字符串(错误信息)，则放到mssages列表中，而原genes.de信息清空
    if (is.character(x = genes.de[[i]])) {
      messages[[i]] <- genes.de[[i]]
      genes.de[[i]] <- NULL
    }
  }


  #(A6) 整理返回结果：新建空数据框
  gde.all <- data.frame()
  

  #(A7) 遍历每个类
  for (i in 1:length(x = idents.all)) {
  	#(B1) 如果是空，则下一个循环
    if (is.null(x = unlist(x = genes.de[i]))) {
      next
    }
    #(B2) 对于临时变量gde
    gde <- genes.de[[i]]
    #(B3) 如果行数>0
    if (nrow(x = gde) > 0) {
      #(C1) 如果是 roc 检验，按myAUC列取子集
      if (test.use == "roc") {
        gde <- subset(
          x = gde,
          subset = (myAUC > return.thresh | myAUC < (1 - return.thresh))
        )
      #(C2) 如果node为空，或者 是 bimod 或 t 检验： 默认wilcox检验也在这个范围
      } else if (is.null(x = node) || test.use %in% c('bimod', 't')) {
      	# 按照(p值 升序，第二列FC降序)排序
        gde <- gde[order(gde$p_val, -gde[, 2]), ]
        # 按p值取子集
        gde <- subset(x = gde, subset = p_val < return.thresh)
      }

      #(C3) 如果剩余行数>0
      if (nrow(x = gde) > 0) {
      	# 新增一列 cluster
        gde$cluster <- idents.all[i]
        # 新增一列 gene
        gde$gene <- rownames(x = gde)
      }

      #(C4) 如果行数>0，合并到 gde.all 数据框中（可以和上一个if合并）
      if (nrow(x = gde) > 0) {
        gde.all <- rbind(gde.all, gde)
      }
    }
  }


  #(A8) 如果定义只要正marker，且行数>0
  if ((only.pos) && nrow(x = gde.all) > 0) {
  	# 按照第二列>0子集，并返回
    return(subset(x = gde.all, subset = gde.all[, 2] > 0))
  }

  #(A9) 取添加的基因列，唯一化，作为行名: 有什么意义？万一两个类都有这个基因怎么办？
  rownames(x = gde.all) <- make.unique(names = as.character(x = gde.all$gene))

  #(A10) 如果行数==0，立刻警告
  if (nrow(x = gde.all) == 0) {
    warning("No DE genes identified", call. = FALSE, immediate. = TRUE)
  }

  #(A11) 如果有(错误)消息：则该类没有比较
  if (length(x = messages) > 0) {
    warning("The following tests were not performed: ", call. = FALSE, immediate. = TRUE)
    for (i in 1:length(x = messages)) {
      if (!is.null(x = messages[[i]])) {
        warning("When testing ", idents.all[i], " versus all:\n\t", messages[[i]], call. = FALSE, immediate. = TRUE)
      }
    }
  }

  #(A12) 如果node非空，求子集: 一般跳过
  if (!is.null(x = node)) {
    gde.all$cluster <- MapVals(
      vec = gde.all$cluster,
      from = new.nodes,
      to = orig.nodes
    )
  }

  #(A13) 返回数据框
  return(gde.all)
}









() FeaturePlot
$ find . | grep -P "(R|c)$" | grep -v "notes"  | xargs grep -in "FeaturePlot" --color=auto 2>/dev/null
./seurat-4.1.0/R/visualization.R:971:FeaturePlot <- function(


#' Visualize 'features' on a dimensional reduction plot
#'
#' Colors single cells on a dimensional reduction plot according to a 'feature'
#' (i.e. gene expression, PC scores, number of genes detected, etc.)
#'
#' @inheritParams DimPlot
#' @param order Boolean determining whether to plot cells in order of expression. Can be useful if
#' cells expressing given feature are getting buried.
#' @param features Vector of features to plot. Features can come from:
#' \itemize{
#'     \item An \code{Assay} feature (e.g. a gene name - "MS4A1")
#'     \item A column name from meta.data (e.g. mitochondrial percentage - "percent.mito")
#'     \item A column name from a \code{DimReduc} object corresponding to the cell embedding values
#'     (e.g. the PC 1 scores - "PC_1")
#' }
#' @param cols The two colors to form the gradient over. Provide as string vector with
#' the first color corresponding to low values, the second to high. Also accepts a Brewer
#' color scale or vector of colors. Note: this will bin the data into number of colors provided.
#' When blend is \code{TRUE}, takes anywhere from 1-3 colors:
#' \describe{
#'   \item{1 color:}{Treated as color for double-negatives, will use default colors 2 and 3 for per-feature expression}
#'   \item{2 colors:}{Treated as colors for per-feature expression, will use default color 1 for double-negatives}
#'   \item{3+ colors:}{First color used for double-negatives, colors 2 and 3 used for per-feature expression, all others ignored}
#' }
#' @param min.cutoff,max.cutoff Vector of minimum and maximum cutoff values for each feature,
#'  may specify quantile in the form of 'q##' where '##' is the quantile (eg, 'q1', 'q10')
#' @param split.by A factor in object metadata to split the feature plot by, pass 'ident'
#'  to split by cell identity'; similar to the old \code{FeatureHeatmap}
#' @param keep.scale How to handle the color scale across multiple plots. Options are:
#' \itemize{
#'   \item{"feature" (default; by row/feature scaling):}{ The plots for each individual feature are scaled to the maximum expression of the feature across the conditions provided to 'split.by'.}
#'   \item{"all" (universal scaling):}{ The plots for all features and conditions are scaled to the maximum expression value for the feature with the highest overall expression.}
#'   \item{NULL (no scaling):}{ Each individual plot is scaled to the maximum expression value of the feature in the condition provided to 'split.by'. Be aware setting NULL will result in color scales that are not comparable between plots.}
#' }
#' @param slot Which slot to pull expression data from?
#' @param blend Scale and blend expression values to visualize coexpression of two features
#' @param blend.threshold The color cutoff from weak signal to strong signal; ranges from 0 to 1.
#' @param ncol Number of columns to combine multiple feature plots to, ignored if \code{split.by} is not \code{NULL}
#' @param coord.fixed Plot cartesian coordinates with fixed aspect ratio
#' @param by.col If splitting by a factor, plot the splits per column with the features as rows; ignored if \code{blend = TRUE}
#' @param sort.cell Redundant with \code{order}. This argument is being
#' deprecated. Please use \code{order} instead.
#' @param interactive Launch an interactive \code{\link[Seurat:IFeaturePlot]{FeaturePlot}}
#' @param combine Combine plots into a single \code{\link[patchwork]{patchwork}ed}
#' ggplot object. If \code{FALSE}, return a list of ggplot objects
#'
#' @return A \code{\link[patchwork]{patchwork}ed} ggplot object if
#' \code{combine = TRUE}; otherwise, a list of ggplot objects
#'
#' @importFrom grDevices rgb
#' @importFrom patchwork wrap_plots
#' @importFrom cowplot theme_cowplot
#' @importFrom RColorBrewer brewer.pal.info
#' @importFrom ggplot2 labs scale_x_continuous scale_y_continuous theme element_rect
#' dup_axis guides element_blank element_text margin scale_color_brewer scale_color_gradientn
#' scale_color_manual coord_fixed ggtitle
#'
#' @export
#' @concept visualization
#'
#' @note For the old \code{do.hover} and \code{do.identify} functionality, please see
#' \code{HoverLocator} and \code{CellSelector}, respectively.
#'
#' @aliases FeatureHeatmap
#' @seealso \code{\link{DimPlot}} \code{\link{HoverLocator}}
#' \code{\link{CellSelector}}
#'
#' @examples
#' data("pbmc_small")
#' FeaturePlot(object = pbmc_small, features = 'PC_1')
#'
FeaturePlot <- function(
  object,
  features,
  dims = c(1, 2),
  cells = NULL,
  cols = if (blend) {
    c('lightgrey', '#ff0000', '#00ff00')
  } else {
    c('lightgrey', 'blue')
  },
  pt.size = NULL,
  order = FALSE,
  min.cutoff = NA,
  max.cutoff = NA,
  reduction = NULL,
  split.by = NULL,
  keep.scale = "feature",
  shape.by = NULL,
  slot = 'data',
  blend = FALSE,
  blend.threshold = 0.5,
  label = FALSE,
  label.size = 4,
  label.color = "black",
  repel = FALSE,
  ncol = NULL,
  coord.fixed = FALSE,
  by.col = TRUE,
  sort.cell = NULL,
  interactive = FALSE,
  combine = TRUE,
  raster = NULL,
  raster.dpi = c(512, 512)
) {
  #(A1) 参数 sort.cell 作废，推荐使用 order
  # TODO: deprecate fully on 3.2.0
  if (!is.null(x = sort.cell)) {
    warning(
      "The sort.cell parameter is being deprecated. Please use the order ",
      "parameter instead for equivalent functionality.",
      call. = FALSE,
      immediate. = TRUE
    )
    #如果使用老参数 sort.cell 是T，则覆盖掉 order 参数
    if (isTRUE(x = sort.cell)) {
      order <- sort.cell
    }
  }

  #(A2) 如果是交互模式，则调用 IFeaturePlot 函数并返回
  if (interactive) {
    return(IFeaturePlot(
      object = object,
      feature = features[1],
      dims = dims,
      reduction = reduction,
      slot = slot
    ))
  }

  #(A3) 检查缩放设置 合法性: 非空，且不在 c("feature", "all") 中，则报错
  # Check keep.scale param for valid entries
  if (!(is.null(x = keep.scale)) && !(keep.scale %in% c("feature", "all"))) {
    stop("`keep.scale` must be set to either `feature`, `all`, or NULL")
  }


  #(A4) 设置主题，删除右边的Y轴的竖线、刻度、刻度文字
  # Set a theme to remove right-hand Y axis lines
  # Also sets right-hand Y axis text label formatting
  no.right <- theme(
    axis.line.y.right = element_blank(),
    axis.ticks.y.right = element_blank(),
    axis.text.y.right = element_blank(),
    axis.title.y.right = element_text(
      face = "bold",
      size = 14,
      margin = margin(r = 7)
    )
  )


  #(A5) 获取降维坐标，默认是：> DefaultDimReduc(pbmc) #[1] "umap"
  # Get the DimReduc to use
  reduction <- reduction %||% DefaultDimReduc(object = object)

  #(A6) 如果dims 长度不是2，或者不是数字：报错
  if (length(x = dims) != 2 || !is.numeric(x = dims)) {
    stop("'dims' must be a two-length integer vector")
  }

  #(A7) 如果设置了 bend=T，且 feature 不是2个，则报错
  # Figure out blending stuff
  if (blend && length(x = features) != 2) {
    stop("Blending feature plots only works with two features")
  }

  #(A8) 如果 bend==T，设置颜色
  # Set color scheme for blended FeaturePlots
  if (blend) {
  	# 获取函数的参数列表的默认值
    default.colors <- eval(expr = formals(fun = FeaturePlot)$cols)
    # 获取用户设置的值 cols
    cols <- switch(
      # 长度 分情况给提示
      EXPR = as.character(x = length(x = cols)),
      '0' = {
        warning(
          "No colors provided, using default colors",
          call. = FALSE,
          immediate. = TRUE
        )
        default.colors #只有0个用默认值
      },
      '1' = {
        warning(
          "Only one color provided, assuming specified is double-negative and augmenting with default colors",
          call. = FALSE,
          immediate. = TRUE
        )
        c(cols, default.colors[2:3]) #只有一个颜色，使用：该颜色 + 默认的后两种颜色
      },
      '2' = {
        warning(
          "Only two colors provided, assuming specified are for features and agumenting with '",
          default.colors[1],
          "' for double-negatives",
          call. = FALSE,
          immediate. = TRUE
        )
        c(default.colors[1], cols) #只有2个颜色，使用默认的第一个颜色 + 这2个颜色
      },
      '3' = cols, #3个直接用

      {
        warning(
          "More than three colors provided, using only first three",
          call. = FALSE,
          immediate. = TRUE
        )
        cols[1:3] # 否则，就是超过3个颜色，只使用前3个
      }
    )
  }

  #(A9) 如果blend ==T且颜色不是3个，则报错
  if (blend && length(x = cols) != 3) {
    stop("Blending feature plots only works with three colors; first one for negative cells")
  }

  #(A10) 降维的名字
  # dims默认c(1,2), 
  # > Key(pbmc[["umap"]]) #[1] "UMAP_"
  # 合并后就是 UMAP_1
  # Name the reductions
  dims <- paste0(Key(object = object[[reduction]]), dims)

  #(A11) 细胞名，默认是所有（列名）
  cells <- cells %||% colnames(x = object)


  #(A12) 获取用于绘图的数据: 这是该包中我最喜欢的函数之一，可以节省很多代码，且速度足够快。
  # Get plotting data
  	#> FetchData(pbmc, vars = c("UMAP_1", "UMAP_2", "ident", "CD79A", "CD8A"), slot="data") |> head(2)
	#> FetchData(pbmc, vars = c("UMAP_1", "UMAP_2", "ident", "CD79A", "CD8A")) |> head(2)
	#                   UMAP_1   UMAP_2 ident    CD79A     CD8A
	#AAACATACAACCAC-1 2.864640   4.0769     2 0.000000 1.635873
	#AAACATTGAGCTAC-1 5.019568 -12.4723     3 1.962726 0.000000
  data <- FetchData( #该函数见 详解9
    object = object,
    vars = c(dims, 'ident', features),
    cells = cells,
    slot = slot
  )


  #(A13) 检查是否存在：特征、降维
  # Check presence of features/dimensions
  # 如果仅检测到3列，则报错
  if (ncol(x = data) < 4) {
    stop(
      "None of the requested features were found: ",
      paste(features, collapse = ', '),
      " in slot ",
      slot,
      call. = FALSE
    )
  #如果不是所有的dims都存在 新获取数据的列名中，则报错
  } else if (!all(dims %in% colnames(x = data))) {
    stop("The dimensions requested were not found", call. = FALSE)
  }

  #(A14) 特征是：列名从第4及之后的。其他没获取到。
  features <- colnames(x = data)[4:ncol(x = data)]

  #(A15) 确定阈值
  # Determine cutoffs
  min.cutoff <- mapply( #mapply 参数1是函数，后面是该函数的 参数列表
    FUN = function(cutoff, feature) {
      return(ifelse(
        test = is.na(x = cutoff),
        yes = min(data[, feature]), #如果cutoff==NA(默认):则返回基因列中的最值
        no = cutoff                #设置过，则使用设置的值
      ))
    },
    cutoff = min.cutoff,
    feature = features #每一个参数列给一个返回值
  )

  max.cutoff <- mapply(
    FUN = function(cutoff, feature) {
      return(ifelse(
        test = is.na(x = cutoff),
        yes = max(data[, feature]),
        no = cutoff
      ))
    },
    cutoff = max.cutoff,
    feature = features
  )

  #(A16) 检查长度: 要求三个长度一致。否则报错
  check.lengths <- unique(x = vapply(
    X = list(features, min.cutoff, max.cutoff),
    FUN = length,
    FUN.VALUE = numeric(length = 1) #vapply可以定义每个返回值的类型：这里是长度为1的数字
  ))
  if (length(x = check.lengths) != 1) {
    stop("There must be the same number of minimum and maximum cuttoffs as there are features")
  }

  #(A17) 获取颜色长度：如果只提供一个颜色，则作为RColorBrewer包的预制配色方案下标
  brewer.gran <- ifelse(
    test = length(x = cols) == 1,
    yes = brewer.pal.info[cols, ]$maxcolors,
    no = length(x = cols)
  )

  #(A18) 应用以上cutoffs值
  # Apply cutoffs
  data[, 4:ncol(x = data)] <- sapply(
    X = 4:ncol(x = data), #遍历值是列下标，从第4列开始
    FUN = function(index) { #对某一列
      data.feature <- as.vector(x = data[, index]) #取该列数据
      #(B1)第一个阈值，下标是 index-3; 如果用户没有设置阈值，则使用的是该列的最值
      min.use <- SetQuantile(cutoff = min.cutoff[index - 3], data.feature)
      max.use <- SetQuantile(cutoff = max.cutoff[index - 3], data.feature)
      # 小于or大于最值的设置为最值
      data.feature[data.feature < min.use] <- min.use
      data.feature[data.feature > max.use] <- max.use
      #(B2) 如果颜色长度是2，则返回该列数据
      if (brewer.gran == 2) {
        return(data.feature)
      }
      #(B3) 超过2个颜色的才会继续，默认2个颜色。
      # 啥时候三个颜色呢？blend=T时

      #(B4)如果所有数据都是0，则截断值是0
      data.cut <- if (all(data.feature == 0)) {
        0
      }
      # 否则，把该列数据按照 颜色个数切分，并转为因子
      else {
        as.numeric(x = as.factor(x = cut(
          x = as.numeric(x = data.feature),
          breaks = brewer.gran
        )))
      }
      # 返回截断分类后的因子
      return(data.cut)
    }
  )

  #(A19) 给第4列及以后的列名，是否有必要？不是已经有了？
  #经第3(10)验证，sapply返回值确实没有行列名
  colnames(x = data)[4:ncol(x = data)] <- features
  rownames(x = data) <- cells



  #(A20) 处理 split.by 参数，用来分面
  # Figure out splits (FeatureHeatmap)
  #  (B1)如果没有 split.by 参数，则随机设置
  data$split <- if (is.null(x = split.by)) {
    RandomName()
  } else {
  	# (B2) 如果设置了 split.by
    switch(
      EXPR = split.by,
      # 如果 split.by=="ident"，则返回所有细胞的idents，按细胞取子集
      ident = Idents(object = object)[cells, drop = TRUE],
      # 否则，取 split.by 列，按细胞取子集
      # > pbmc[["nFeature_RNA", drop = TRUE]] |> head(3)
	  #AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 
	  #             779             1352             1129 
      object[[split.by, drop = TRUE]][cells, drop = TRUE]
    )
  }
  # split 变量如果不是factor，转为factor
  if (!is.factor(x = data$split)) {
    data$split <- factor(x = data$split)
  }


  #(A21) 形状参数
  # Set shaping variable
  if (!is.null(x = shape.by)) {
    data[, shape.by] <- object[[shape.by, drop = TRUE]]
  }

  #(A22) 通常不适用blend，所以 length = features数 * split的因子个数
  # Make list of plots
  plots <- vector(
    mode = "list",
    length = ifelse(
      test = blend,
      yes = 4, #2个是x和y，剩下2列是基因
      no = length(x = features) * length(x = levels(x = data$split))
    )
  )

  #(A23) 应用常规限制 Apply common limits
  # x和y轴的范围：取UMAP_1列的最小值，向下取整数；该列最大值，向上取整数
  xlims <- c(floor(x = min(data[, dims[1]])), ceiling(x = max(data[, dims[1]])))
  #大佬的代码风格也不一致，和上一行的floor()比，没有写 x=
  ylims <- c(floor(min(data[, dims[2]])), ceiling(x = max(data[, dims[2]]))) 

  #(A24) 设置 blend 的颜色: 跳过
  # Set blended colors
  if (blend) {
    ncol <- 4
    color.matrix <- BlendMatrix(
      two.colors = cols[2:3],
      col.threshold = blend.threshold,
      negative.color = cols[1]
    )

    cols <- cols[2:3]

    colors <- list(
      color.matrix[, 1],
      color.matrix[1, ],
      as.vector(x = color.matrix)
    )
  }


  #(A25) 散点图
  # Make the plots
  # 遍历每个 split。没有给split.by参数也默认有一个临时5字母名字
  for (i in 1:length(x = levels(x = data$split))) {
    #(B1) 需要绘制哪个 split
    # Figure out which split we're working with
    ident <- levels(x = data$split)[i]

    #(B2) 取该ident的子集，保持df形状
    # 注意：这里比较因子时，先转为字符串！
    data.plot <- data[as.character(x = data$split) == ident, , drop = FALSE]

    #(B3) blend：跳过
    # Blend expression values
    if (blend) {
      features <- features[1:2]
      no.expression <- features[colMeans(x = data.plot[, features]) == 0]
      if (length(x = no.expression) != 0) {
        stop(
          "The following features have no value: ",
          paste(no.expression, collapse = ', '),
          call. = FALSE
        )
      }
      data.plot <- cbind(data.plot[, c(dims, 'ident')], BlendExpression(data = data.plot[, features[1:2]]))
      features <- colnames(x = data.plot)[4:ncol(x = data.plot)]
    }

    #(B4) 遍历 feature
    # Make per-feature plots
    for (j in 1:length(x = features)) {
      #(C1) 获取1个基因名
      feature <- features[j]

      #(C2) blend 跳过； False 时，cols.use=NULL
      # Get blended colors
      if (blend) {
        cols.use <- as.numeric(x = as.character(x = data.plot[, feature])) + 1
        cols.use <- colors[[j]][sort(x = unique(x = cols.use))]
      } else {
        cols.use <- NULL
      }

      #(C3) 取单个数据：x和y，ident列，一个基因，形状(默认 NULL)
      data.single <- data.plot[, c(dims, 'ident', feature, shape.by)]
      # Make the plot
      plot <- SingleDimPlot(  #该函数见 解析24-2(2)
        data = data.single,
        dims = dims,
        col.by = feature,
        order = order,
        pt.size = pt.size,
        cols = cols.use,
        shape.by = shape.by,
        label = FALSE,
        raster = raster,
        raster.dpi = raster.dpi
      ) +
        #限定x和y绘制范围
        scale_x_continuous(limits = xlims) +
        scale_y_continuous(limits = ylims) +
        theme_cowplot() +
        CenterTitle() #标题居中
        # theme(plot.title = element_text(hjust = 0.5))
      # Add labels

      #(C4) 如果label=T，标记上 ident 列信息
      if (label) {
        plot <- LabelClusters( #见 解析24-2(4)
          plot = plot,
          id = 'ident',
          repel = repel,
          size = label.size,
          color = label.color
        )
      }

      #(C5) 美化：如果split因子超过1个。一般指定 split.by 参数，就走这一步
      # Make FeatureHeatmaps look nice(ish)
      if (length(x = levels(x = data$split)) > 1) {
      	# (D1)分面标题：没有边框，字体黑色
        plot <- plot + theme(panel.border = element_rect(fill = NA, colour = 'black'))
        # Add title 添加标题：只有ident==1的图添加title
        plot <- plot + if (i == 1) {
          labs(title = feature)
        } else {
          labs(title = NULL)
        }

        # (D2)添加第二个轴，如果j是本ident的最后一个 feature，且 无blend
        # Add second axis
        if (j == length(x = features) && !blend) {
          suppressMessages(
            expr = plot <- plot +
              # 添加第二条y坐标轴
              scale_y_continuous(
                sec.axis = dup_axis(name = ident),
                limits = ylims
              ) +
              no.right #但是又不显示右侧坐标轴。这个又啥用？ //todo
          )
        }

        # (D3)如果不是第一个 feature，去掉左侧的y轴
        # Remove left Y axis
        if (j != 1) {
          plot <- plot + theme(
            axis.line.y = element_blank(),
            axis.ticks.y = element_blank(),
            axis.text.y = element_blank(),
            axis.title.y.left = element_blank()
          )
        }

        # (D4)如果不是最后一个 ident，去掉x轴
        # Remove bottom X axis
        if (i != length(x = levels(x = data$split))) {
          plot <- plot + theme(
            axis.line.x = element_blank(),
            axis.ticks.x = element_blank(),
            axis.text.x = element_blank(),
            axis.title.x = element_blank()
          )
        }
      #(C6) 如果没有 split.by 参数，默认不分面
      } else {
      	# 加标题：feature
        plot <- plot + labs(title = feature)
      }

      #(C7) 添加颜色映射，非 blend 模式(大多数走这一步)
      # Add colors scale for normal FeaturePlots
      if (!blend) {
      	# 这一行多余? //todo
        plot <- plot + guides(color = NULL)
        # (D1)渐变色：记录参数
        cols.grad <- cols

        # (D2)如果1个，则使用调色板 //todo 我尝试调用色板，失败。只是当作单一颜色使用。
        if (length(x = cols) == 1) {
          plot <- plot + scale_color_brewer(palette = cols)
        # (D3)默认给出2个颜色
        } else if (length(x = cols) > 1) {
          # 取基因列的uniq值
          unique.feature.exp <- unique(data.plot[, feature])
          # 如果只有1个值，可能表达相同，或者都是0
          if (length(unique.feature.exp) == 1) {
            warning("All cells have the same value (", unique.feature.exp, ") of ", feature, ".")
            # 如果都相等，且等于0，则渐变色是第1个：背景色
            if (unique.feature.exp == 0) {
              cols.grad <- cols[1]
            # 都相等且不等于0，则是颜色数组，有啥意义？ //todo
            } else{
              cols.grad <- cols
            }
          }
          # 使用渐变色，默认2色走这一步
          plot <- suppressMessages(
            expr = plot + scale_color_gradientn(
              colors = cols.grad, #来自C7-D1，然后又被覆盖一次
              guide = "colorbar"
            )
          )
        }
      }


      # (C8) 如果 keep.scale 非空，且为 feature，且 非blend：默认走这一步
      if (!(is.null(x = keep.scale)) && keep.scale == "feature" && !blend) {
      	# 获取该列最值
        max.feature.value <- max(data[, feature])
        min.feature.value <- min(data[, feature])
        # 颜色又被覆盖一次，限定颜色的范围为该基因的2个最值之间
        plot <- suppressMessages(plot & 
        	scale_color_gradientn(
        		colors = cols,
        		limits = c(min.feature.value, max.feature.value)))
      }

      #(C9) 是否固定x和y的比例，默认F
      # Add coord_fixed
      if (coord.fixed) {
        plot <- plot + coord_fixed()
      }

      #(C10) 魔术方法：作者认为不这么写有时候不行！不知道为啥
      # I'm not sure why, but sometimes the damn thing fails without this
      # Thanks ggplot2
      plot <- plot

      #(C11) 把小图放到list中。下标: 外循环i已经满的行数 * 基因数 + 当前基因数
      # i 是 外循环 ident 分面，
      #    j是内循环 基因
      # Place the plot
      plots[[(length(x = features) * (i - 1)) + j]] <- plot
    }

  }





  #(A26) blend 模式的颜色key: 跳过
  # Add blended color key
  if (blend) {
    blend.legend <- BlendMap(color.matrix = color.matrix)
    for (ii in 1:length(x = levels(x = data$split))) {
      suppressMessages(expr = plots <- append(
        x = plots,
        values = list(
          blend.legend +
            scale_y_continuous(
              sec.axis = dup_axis(name = ifelse(
                test = length(x = levels(x = data$split)) > 1,
                yes = levels(x = data$split)[ii],
                no = ''
              )),
              expand = c(0, 0)
            ) +
            labs(
              x = features[1],
              y = features[2],
              title = if (ii == 1) {
                paste('Color threshold:', blend.threshold)
              } else {
                NULL
              }
            ) +
            no.right
        ),
        after = 4 * ii - 1
      ))
    }
  }


  #(A27) 去掉NUL的散点图：啥时候会被去掉？想不到。因为前面 FetchData 已经算一次过滤了。
  # Remove NULL plots
  plots <- Filter(f = Negate(f = is.null), x = plots)

  #(A28) 合并图形
  # Combine the plots
  if (is.null(x = ncol)) { # 如果没有设置 ncol，则默认为2
    ncol <- 2
    if (length(x = features) == 1) { #如果只有1个基因，则 ncol=1
      ncol <- 1
    }
    if (length(x = features) > 6) { # 如果超过6个基因，则 ncol=3
      ncol <- 3
    }
    if (length(x = features) > 9) { #超过9个基因，则 ncol=4
      ncol <- 4
    }
  }

  #(A28) 设置列数
  # 如果没有split，或者是blend模式，ncol不变。
  # 否则：split.by非空时，ncol=features个数
  ncol <- ifelse(
    test = is.null(x = split.by) || blend,
    yes = ncol,
    no = length(x = features)
  )

  #(A29) 如果 blend 模式：legend='none';
  # 非 blend 模式：如果 split.by不为null，则返回'none'； 如果为null，则返回null。
  # 像是一句废话： null和 'none' 的区别是啥？ //todo
  legend <- if (blend) {
    'none'
  } else {
    split.by %iff% 'none'
  }

  #(A30) 如果要合并（默认），仅对 非blend 模式
  # Transpose the FeatureHeatmap matrix (not applicable for blended FeaturePlots)
  if (combine) {
  	# (B1) 按列split(by.col默认T)，且 spli.by非空，非blend模式
    if (by.col && !is.null(x = split.by) && !blend) {

      #(C1) 对list的每个元素使用FUN函数：添加主题，标题，y轴显示第二个轴+范围+不显示右侧坐标轴(这个有意义吗？)
      plots <- lapply(
        X = plots,
        FUN = function(x) {
          return(suppressMessages(
            expr = x +
              theme_cowplot() +
              ggtitle("") +
              scale_y_continuous(sec.axis = dup_axis(name = ""), limits = ylims) +
              no.right
          ))
        }
      )

      #(C2) 获取 split 的个数
      nsplits <- length(x = levels(x = data$split))

      #
      # 经过代码前后印证，默认，这些plot的编号是按列by.col=T
      # 比如2个基因（行），split.by分成3列，共6个小图。编号分别是第一列12，第二列34，第三列56
      # 1 3 5
      # 2 4 6
      # 所以，5，6是最后一列，添加基因名
      # 1 3 5 是第一行，添加分面名
      

      #(C3) 遍历最后一列: 上一列结尾+1，到这一列结尾。最后一列添加基因名字
      ## 2基因(行) * 3列
	  #	for (i in (2 * (3 - 1) + 1):(2 * 3)) {
	  #	  print(i)
	  #	}
	  # [1] 5
	  # [1] 6
      idx <- 1 #idx=1，基因(features)的下标
      for (i in (length(x = features) * (nsplits - 1) + 1):(length(x = features) * nsplits)) {
        plots[[i]] <- suppressMessages(
          # 这个有啥意义？把C1的lapply换for循环又搞了一遍。测试：之前C1去掉，则左下角图又有title，不美观
          # 之前是去掉名字，这里是加上基因名，最后一列。
          expr = plots[[i]] +
            scale_y_continuous(
              sec.axis = dup_axis(name = features[[idx]]), #加上基因名字
              # 为啥用[[]]而不是[]?就是写错了，数组取元素[]，向量取元素[[]]。虽然[[]]用于数组取元素结果没区别。
              limits = ylims
            ) +
            no.right
        )
        idx <- idx + 1
      }


      #(C4) 遍历小图：对于第一行的图，添加分面名
	  # > 1:6 %% 2 #对于6个图，2个基因（一行一个基因）：3列
	  #[1] 1 0 1 0 1 0
	  #> 1:6 %% 2 == 1
	  #[1]  TRUE FALSE  TRUE FALSE  TRUE FALSE
	  #> which(1:6 %% 2 == 1)
	  #[1] 1 3 5
      idx <- 1
      for (i in which(x = 1:length(x = plots) %% length(x = features) == 1)) {
        plots[[i]] <- plots[[i]] +
          ggtitle(levels(x = data$split)[[idx]]) +
          theme(plot.title = element_text(hjust = 0.5)) #标题居中
        idx <- idx + 1
      }


      #(C5) 如果只有一个基因(仅一行)，修正C4中 features 只有一个，导致无法进入循环的bug:
      # > which(1:6 %% 1 == 1)
      # integer(0)
      idx <- 1
      if (length(x = features) == 1) {
      	# 对于每个小图
        for (i in 1:length(x = plots)) {
          plots[[i]] <- plots[[i]] +
            ggtitle(levels(x = data$split)[[idx]]) +
            theme(plot.title = element_text(hjust = 0.5))
          idx <- idx + 1
        }
        # 一列，行数按 split：我的测试表明，是1行 2列 //todo
        # FeaturePlot2(object = pbmc, features = c( "CD3G"), 
        #       split.by = "gene_high", pt.size = 1)
        ncol <- 1
        nrow <- nsplits
      } else {
      #(C6) 有多个基因时，行数等于：如果split.by定义了，则使用后面的 split 因子个数。
      # 如果 split.by==NULL(默认)，则 nrow= NULL
        nrow <- split.by %iff% length(x = levels(x = data$split))
      }
      #todo: split不应该是列吗？ 一行是一个基因


      #(C7) 调整小图的顺序，代码压缩度很高，详见本文3(27)
      plots <- plots[c(do.call(
        what = rbind,
        args = split(
        		x = 1:length(x = plots), 
        		f = ceiling(x = seq_along(along.with = 1:length(x = plots)) / length(x = features))
        	)
      ))]
      
      #(C8) 果然，这里行列又颠倒回去了....：行=列，列=行!?
      # 不知道作者怎么想的，反正我感觉前面反了，这里再返回去。 为什么要这样？
      # 分析： 绘制plots时，两层循环，外是split，内是gene。内循环变化快，外循环变化慢。
      # Set ncol to number of splits (nrow) and nrow to number of features (ncol)
      plots <- wrap_plots(plots, ncol = nrow, nrow = ncol)

      #(C9) 是否加图例：如果 legend 非空，且 legend =='none'
      #往前找，legend不是函数参数，取决于 split.by 参数: 通常blend=F，legend=split.by %iff% 'none'
	      # 如果 split.by=null, legend=null；
	      # 如果 split.by非空， legend='none'
      if (!is.null(x = legend) && legend == 'none') {
        plots <- plots & NoLegend() #去掉图例
        # 什么情况下有图例呢？ split.by=NULL时
      }
    #(B2) 其他情况
    # * 按列split(by.col默认T)，如果 split.by 为空。
    # * 或blend模式。
    # * 或 by.col=F。
    } else {
      plots <- wrap_plots(plots, 
      	ncol = ncol, 
      	nrow = split.by %iff% length(x = levels(x = data$split)))
    }

    #(B3) 又来一遍 B1-C9
    if (!is.null(x = legend) && legend == 'none') {
      plots <- plots & NoLegend()
    }

    #(B4) 如果 keep.scale非空，且值为 all（默认不是这个值），非blend模式
    if (!(is.null(x = keep.scale)) && keep.scale == "all" && !blend) {
      # 找到全局最值，并应用到颜色中
      max.feature.value <- max(data[, features])
      min.feature.value <- min(data[, features])
      plots <- suppressMessages(plots & 
      	scale_color_gradientn(
      		colors = cols, 
      		limits = c(min.feature.value, max.feature.value))
      	)
    }
  }

  #返回。默认是 拼合后的 ggplot2对象
  return(plots)
}






() SetQuantile
./seurat-4.1.0/R/utilities.R:1527:SetQuantile <- function(cutoff, data) {

#' Find the Quantile of Data
寻找分位数，函数名不应该是 get 吗？
#'
#' Converts a quantile in character form to a number regarding some data.
#' String form for a quantile is represented as a number prefixed with
#' \dQuote{q}; for example, 10th quantile is \dQuote{q10} while 2nd quantile is
#' \dQuote{q2}. Will only take a quantile of non-zero data values
把分位数字符转为数字，比如 q10 表示10th分位数；q2表示2nd分位数。
只对非0数字计算分位数。
#'
#' @param cutoff The cutoff to turn into a quantile
#' @param data The data to turn find the quantile of
#'
#' @return The numerical representation of the quantile
#'
#' @importFrom stats quantile
#'
#' @export
#' @concept utilities
#'
#' @examples
#' set.seed(42)
#' SetQuantile('q10', sample(1:100, 10))
#'
SetQuantile <- function(cutoff, data) {
  #cutoff的形式: 开头是q，然后是1-2位数字0-9。也就是说 q100是不行的。
  if (grepl(pattern = '^q[0-9]{1,2}$', x = as.character(x = cutoff), perl = TRUE)) {
    # 替换q为空，转为数值，然后除以100：q10-> 10 -> 10/100=0.1
    this.quantile <- as.numeric(x = sub(
      pattern = 'q',
      replacement = '',
      x = as.character(x = cutoff)
    )) / 100

    # 对数据解list
    data <- unlist(x = data)
    # 只取正数
    data <- data[data > 0]
    
    #对数据求分位数 
    cutoff <- quantile(x = data, probs = this.quantile)
  }
  # 否则，直接返回该 cutoff 强转数值的值
  return(as.numeric(x = cutoff))
}



